Data prep

suppressWarnings(library(MendelianRandomization))
suppressWarnings(library(digest))
suppressWarnings(library(tidyr))                         
suppressWarnings(library(MRPRESSO))
suppressWarnings(library(rowr))
suppressWarnings(library(ggplot2))
suppressWarnings(library(dplyr))

setwd('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants')
diabetes_all<-read.csv('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants/Diabetes/Newrf_diabetes_instruments.csv', header = TRUE, stringsAsFactors = FALSE)
#diabetes_all<-diabetes_all[diabetes_all$F.statistic>10,] Use command to remove weak instruments


breast<-readRDS('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants/Breast Cancer/breast.rds')
prostate<-readRDS('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants/Prostate Cancer/prostate.rds')
breast= separate(data=breast, col=phase3_1kg_id, into=c('SNP','rest'), sep= "\\:") #Create SNP column in breast data

#Merge diabetes and breast data
index_b<-match(breast$SNP,diabetes_all$SNP..rs.id.)
diabetes_b<-cbind(diabetes_all[index_b,],breast)

#ALign effect alleles in 2 datasets
data_in_ea = diabetes_b$Effect.Allele
data_in_nea = diabetes_b$Other.Allele

ea_new = diabetes_b$a1
nea_new = diabetes_b$a0


stay=as.numeric(data_in_ea == ea_new)  #stay = 1 if the same
swap=as.numeric(data_in_nea == ea_new) 

#checks
table(stay)
## stay
##  0  1 
## 57 54
table(swap)
## swap
##  0  1 
## 54 57
stay*swap == 0
##   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [57] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [71] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [85] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [99] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
swap_vec = rep(1, length(swap))
swap_vec = (-2 * swap) + swap_vec

#align the beta vector in the new data
diabetes_b$bcac_onco_icogs_gwas_beta =swap_vec*as.numeric(diabetes_b$bcac_onco_icogs_gwas_beta)
diabetes_b$Effect = swap_vec * as.numeric(diabetes_b$Effect)
####################################################################################

#Merge diabetes and prostate data

index_p<-match(prostate$SNP,diabetes_all$SNP..rs.id.)
diabetes_p<-cbind(diabetes_all[index_p,],prostate)

colnames(diabetes_p)[23]<-'Effect_p'
colnames(diabetes_p)[24]<-'StdErr_p'

#ALign effect alleles in 2 datasets
data_in_ea = diabetes_p$Effect.Allele
data_in_nea = diabetes_p$Other.Allele

ea_new = diabetes_p$Allele1
nea_new = diabetes_p$Allele2


stay=as.numeric(data_in_ea == ea_new)  #stay = 1 if the same
swap=as.numeric(data_in_nea == ea_new) 

#checks
table(stay)
## stay
##  0 
## 88
table(swap)
## swap
##  0 
## 88
stay*swap == 0
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [57] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [71] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [85] TRUE TRUE TRUE TRUE
swap_vec = rep(1, length(swap))
swap_vec = (-2 * swap) + swap_vec

#align the beta vector in the new data
diabetes_p$Effect_p =swap_vec*as.numeric(diabetes_p$Effect_p)
diabetes_p$Effect = swap_vec * as.numeric(diabetes_p$Effect)
# remove NAs, we believe that they are missing at random
diabetes_b<-na.omit(diabetes_b)
diabetes_p<-na.omit(diabetes_p)

#Get rid of duplicated columns
diabetes_b <- diabetes_b[, !duplicated(colnames(diabetes_b))]
diabetes_p <- diabetes_p[, !duplicated(colnames(diabetes_p))]
#Load data of rs ids of grouped SNPs from PhenoScanner query

BMI=read.csv('C:/Users/marta/Desktop/IC/Thesis/BMI.csv', stringsAsFactors = FALSE, header=FALSE)
Glucose=read.csv('C:/Users/marta/Desktop/IC/Thesis/Glucose.csv',stringsAsFactors = FALSE,header=FALSE)
Insulin=read.csv('C:/Users/marta/Desktop/IC/Thesis/Insulin.csv',stringsAsFactors = FALSE, header=FALSE)
Lipids=read.csv('C:/Users/marta/Desktop/IC/Thesis/Lipids.csv',stringsAsFactors = FALSE, header=FALSE)

MR breast cancer

MRInputObject_breast <- mr_input(bx = diabetes_b$Effect,
                          bxse = diabetes_b$StdErr,
                          by = diabetes_b$bcac_onco_icogs_gwas_beta,
                          byse = diabetes_b$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=diabetes_b$SNP)

MRAllObject_all_breast <- mr_allmethods(MRInputObject_breast, method = "all")

Breast cancer MR effect estimates:

MRAllObject_all_breast
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median   -0.004     0.019  -0.041 0.033   0.835
##            Weighted median    0.011     0.018  -0.024 0.045   0.539
##  Penalized weighted median   -0.004     0.018  -0.039 0.031   0.825
##                                                                    
##                        IVW    0.000     0.017  -0.033 0.033   0.998
##              Penalized IVW   -0.007     0.013  -0.033 0.019   0.581
##                 Robust IVW    0.005     0.023  -0.040 0.050   0.819
##       Penalized robust IVW   -0.008     0.016  -0.040 0.024   0.612
##                                                                    
##                   MR-Egger    0.029     0.035  -0.040 0.098   0.407
##                (intercept)   -0.003     0.003  -0.009 0.003   0.344
##         Penalized MR-Egger    0.003     0.031  -0.058 0.064   0.926
##                (intercept)   -0.001     0.002  -0.006 0.004   0.743
##            Robust MR-Egger    0.046     0.066  -0.083 0.176   0.484
##                (intercept)   -0.004     0.005  -0.013 0.006   0.447
##  Penalized robust MR-Egger   -0.005     0.041  -0.085 0.076   0.908
##                (intercept)    0.000     0.003  -0.006 0.006   0.939
mr_plot(MRInputObject_breast,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast)

MR-PRESSO

#MR-PRESSO breast

input = data.frame(betaY=diabetes_b$bcac_onco_icogs_gwas_beta, seY=diabetes_b$bcac_onco_icogs_gwas_se,beta_diabetes=diabetes_b$Effect, se_diabetes=diabetes_b$StdErr,row.names= NULL)

mrpresso_out = mr_presso(BetaOutcome = "betaY", BetaExposure = c("beta_diabetes"), SdOutcome = "seY", SdExposure="se_diabetes", data=input, OUTLIERtest = TRUE, DISTORTIONtest = TRUE, NbDistribution = 1000)
## Warning in mr_presso(BetaOutcome = "betaY", BetaExposure =
## c("beta_diabetes"), : Outlier test unstable. The significance threshold
## of 0.05 for the outlier test is not achievable with only 1000 to
## compute the null distribution. The current precision is <0.111. Increase
## NbDistribution.
mrpresso_out$`MR-PRESSO results`$`Global Test`$`RSSobs`
## [1] 459.3003
mrpresso_out$`MR-PRESSO results`$`Global Test`$Pvalue
## [1] "<0.001"
mrpresso_out$`Main MR results`
##        Exposure       MR Analysis Causal Estimate         Sd       T-stat
## 1 beta_diabetes               Raw    4.022679e-05 0.01704120  0.002360562
## 2 beta_diabetes Outlier-corrected   -7.200820e-03 0.01400442 -0.514181786
##     P-value
## 1 0.9981208
## 2 0.6082263

MR-PRESSO Distortion test

mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Distortion Coefficient`
## beta_diabetes 
##      100.5586
mrpresso_out$`MR-PRESSO results`$`Distortion Test`$Pvalue
## [1] 0.228

MR-PRESSO Outliers detected

out_indices = mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Outliers Indices`
outliers = diabetes_b$SNP[out_indices]

# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%outliers,yes='red',no='blue'))) + 
  geom_point(size=1.5) +
  labs(title = "", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "outliers"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 10), axis.title.x = element_text(size = 12),
        axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 12),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

Stratified MR analyses

BMI

Scatter plot of grouped SNPs

# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%BMI$V2,yes='red',no='blue'))) + 
  geom_point(size=1.5) +
  labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "BMI SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

index_bmi<-match(diabetes_b$SNP,BMI$V2)
data_b_bmi<-na.omit(diabetes_b[index_bmi,])
MRInputObject_breast_bmi <- mr_input(bx = data_b_bmi$Effect,
                          bxse = data_b_bmi$StdErr,
                          by = data_b_bmi$bcac_onco_icogs_gwas_beta,
                          byse = data_b_bmi$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_bmi$SNP)

MRAllObject_all_breast_bmi<- mr_allmethods(MRInputObject_breast_bmi, method = "all")
MRAllObject_all_breast_bmi
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median   -0.034     0.029  -0.091 0.023   0.242
##            Weighted median   -0.018     0.027  -0.070 0.035   0.508
##  Penalized weighted median   -0.028     0.027  -0.082 0.025   0.299
##                                                                    
##                        IVW   -0.010     0.025  -0.059 0.039   0.696
##              Penalized IVW   -0.031     0.021  -0.072 0.010   0.138
##                 Robust IVW   -0.010     0.031  -0.071 0.050   0.739
##       Penalized robust IVW   -0.031     0.025  -0.080 0.019   0.224
##                                                                    
##                   MR-Egger    0.052     0.055  -0.056 0.161   0.344
##                (intercept)   -0.006     0.004  -0.014 0.003   0.209
##         Penalized MR-Egger    0.064     0.048  -0.031 0.158   0.188
##                (intercept)   -0.006     0.004  -0.014 0.002   0.128
##            Robust MR-Egger    0.058     0.065  -0.070 0.186   0.375
##                (intercept)   -0.006     0.005  -0.016 0.004   0.215
##  Penalized robust MR-Egger    0.066     0.055  -0.042 0.174   0.233
##                (intercept)   -0.006     0.004  -0.014 0.002   0.134
mr_plot(MRInputObject_breast_bmi,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_bmi)

Insulin

Scatter plot of grouped SNPs

# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%Insulin$V2,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "Insulin SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

index<-match(diabetes_b$SNP,Insulin$V2)
data_b_insulin<-na.omit(diabetes_b[index,])
MRInputObject_breast_insulin <- mr_input(bx = data_b_insulin$Effect,
                          bxse = data_b_insulin$StdErr,
                          by = data_b_insulin$bcac_onco_icogs_gwas_beta,
                          byse = data_b_insulin$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_insulin$SNP)

MRAllObject_all_breast_insulin <- mr_allmethods(MRInputObject_breast_insulin, method = "all")
MRAllObject_all_breast_insulin
##                     Method Estimate Std Error 95% CI         P-value
##              Simple median   -0.072     0.069  -0.207  0.063   0.297
##            Weighted median    0.082     0.047  -0.010  0.175   0.081
##  Penalized weighted median    0.132     0.047   0.040  0.224   0.005
##                                                                     
##                        IVW    0.026     0.062  -0.095  0.146   0.674
##              Penalized IVW   -0.043     0.047  -0.135  0.048   0.350
##                 Robust IVW    0.027     0.061  -0.093  0.146   0.663
##       Penalized robust IVW   -0.041     0.064  -0.167  0.085   0.520
##                                                                     
##                   MR-Egger    0.203     0.111  -0.015  0.422   0.068
##                (intercept)   -0.014     0.008  -0.029  0.001   0.064
##         Penalized MR-Egger    0.284     0.076   0.134  0.434   0.000
##                (intercept)   -0.020     0.005  -0.030 -0.010   0.000
##            Robust MR-Egger    0.221     0.104   0.016  0.425   0.034
##                (intercept)   -0.016     0.007  -0.030 -0.002   0.027
##  Penalized robust MR-Egger    0.286     0.080   0.129  0.443   0.000
##                (intercept)   -0.020     0.006  -0.032 -0.008   0.001
mr_plot(MRInputObject_breast_insulin,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_insulin)

Glucose

Scatter plot of grouped SNPs

# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%Glucose$V2,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "Glucose SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

index<-match(diabetes_b$SNP,Glucose$V2)
data_b_glucose<-na.omit(diabetes_b[index,])
MRInputObject_breast_glucose <- mr_input(bx = data_b_glucose$Effect,
                          bxse = data_b_glucose$StdErr,
                          by = data_b_glucose$bcac_onco_icogs_gwas_beta,
                          byse = data_b_glucose$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_glucose$SNP)

MRAllObject_all_breast_glucose <- mr_allmethods(MRInputObject_breast_glucose, method = "all")
MRAllObject_all_breast_glucose
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median   -0.072     0.076  -0.221 0.078   0.346
##            Weighted median   -0.001     0.058  -0.114 0.113   0.989
##  Penalized weighted median   -0.070     0.061  -0.189 0.049   0.248
##                                                                    
##                        IVW   -0.004     0.071  -0.143 0.136   0.960
##              Penalized IVW   -0.078     0.045  -0.167 0.011   0.087
##                 Robust IVW   -0.010     0.080  -0.168 0.147   0.897
##       Penalized robust IVW   -0.084     0.068  -0.218 0.050   0.220
##                                                                    
##                   MR-Egger    0.067     0.135  -0.198 0.331   0.621
##                (intercept)   -0.005     0.008  -0.022 0.011   0.536
##         Penalized MR-Egger    0.025     0.097  -0.165 0.216   0.794
##                (intercept)   -0.003     0.006  -0.016 0.009   0.610
##            Robust MR-Egger    0.063     0.135  -0.203 0.328   0.642
##                (intercept)   -0.005     0.009  -0.023 0.012   0.539
##  Penalized robust MR-Egger    0.009     0.162  -0.309 0.327   0.955
##                (intercept)   -0.003     0.008  -0.019 0.012   0.690
mr_plot(MRInputObject_breast_glucose,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_glucose)

Lipids

index<-match(diabetes_b$SNP,Lipids$V2)
data_b_lipid<-na.omit(diabetes_b[index,])

Scatter plot of grouped SNPs

# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%Lipids$V2,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "Lipid SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

MRInputObject_breast_lipid <- mr_input(bx = data_b_lipid$Effect,
                          bxse = data_b_lipid$StdErr,
                          by = data_b_lipid$bcac_onco_icogs_gwas_beta,
                          byse = data_b_lipid$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_lipid$SNP)

MRAllObject_all_breast_lipid <- mr_allmethods(MRInputObject_breast_lipid, method = "all")
MRAllObject_all_breast_lipid
##                     Method Estimate Std Error 95% CI         P-value
##              Simple median   -0.077     0.089  -0.252  0.099   0.390
##            Weighted median   -0.075     0.066  -0.205  0.055   0.257
##  Penalized weighted median   -0.130     0.065  -0.257 -0.003   0.045
##                                                                     
##                        IVW   -0.011     0.080  -0.168  0.146   0.894
##              Penalized IVW   -0.105     0.049  -0.202 -0.009   0.033
##                 Robust IVW   -0.025     0.109  -0.239  0.189   0.820
##       Penalized robust IVW   -0.130     0.064  -0.256 -0.005   0.042
##                                                                     
##                   MR-Egger    0.059     0.144  -0.223  0.342   0.680
##                (intercept)   -0.005     0.009  -0.023  0.012   0.554
##         Penalized MR-Egger    0.012     0.104  -0.192  0.217   0.905
##                (intercept)   -0.004     0.007  -0.017  0.010   0.595
##            Robust MR-Egger    0.052     0.150  -0.242  0.345   0.731
##                (intercept)   -0.006     0.009  -0.023  0.012   0.541
##  Penalized robust MR-Egger   -0.053     0.220  -0.484  0.379   0.811
##                (intercept)   -0.002     0.007  -0.017  0.012   0.742
mr_plot(MRInputObject_breast_lipid,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_lipid)

No BMI

BMI=read.csv('C:/Users/marta/Desktop/IC/Thesis/Genetic variants/BMI.csv', stringsAsFactors = FALSE, header=FALSE)
index<-match(diabetes_b$SNP,BMI$V1)
data_b_bmi<-cbind(diabetes_b,BMI[index,])
data_b_bmi$V2[is.na(data_b_bmi$V2)] <- 'NO BMI'
data_nobmi<-data_b_bmi[data_b_bmi$V2=='NO BMI',]

MRInputObject_breast_nobmi <- mr_input(bx = data_nobmi$Effect,
                          bxse = data_nobmi$StdErr,
                          by = data_nobmi$bcac_onco_icogs_gwas_beta,
                          byse = data_nobmi$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_nobmi$SNP)

MRAllObject_all_breast_nobmi <- mr_allmethods(MRInputObject_breast_nobmi, method = "all")
MRAllObject_all_breast_nobmi
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median    0.036     0.030  -0.022 0.095   0.221
##            Weighted median   -0.005     0.027  -0.059 0.049   0.856
##  Penalized weighted median   -0.004     0.028  -0.058 0.050   0.871
##                                                                    
##                        IVW   -0.012     0.022  -0.056 0.031   0.583
##              Penalized IVW    0.002     0.021  -0.038 0.043   0.906
##                 Robust IVW   -0.010     0.022  -0.052 0.033   0.649
##       Penalized robust IVW    0.003     0.020  -0.036 0.041   0.898
##                                                                    
##                   MR-Egger   -0.040     0.048  -0.134 0.055   0.408
##                (intercept)    0.002     0.003  -0.004 0.009   0.517
##         Penalized MR-Egger   -0.043     0.046  -0.133 0.048   0.358
##                (intercept)    0.003     0.003  -0.004 0.009   0.374
##            Robust MR-Egger   -0.047     0.041  -0.127 0.034   0.254
##                (intercept)    0.003     0.004  -0.004 0.011   0.400
##  Penalized robust MR-Egger   -0.049     0.039  -0.126 0.029   0.217
##                (intercept)    0.004     0.003  -0.003 0.010   0.296
mr_plot(MRInputObject_breast_nobmi,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_nobmi)

Prostate cancer MR effect estimates:

MRInputObject_prostate <- mr_input(bx = diabetes_p$Effect,
                          bxse = diabetes_p$StdErr,
                          by = diabetes_p$Effect_p,
                          byse = diabetes_p$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=diabetes_p$SNP)

MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate, method = "all")
MRAllObject_all_prostate 
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median    0.011     0.024  -0.035 0.058   0.632
##            Weighted median    0.001     0.022  -0.042 0.045   0.947
##  Penalized weighted median   -0.005     0.022  -0.049 0.039   0.819
##                                                                    
##                        IVW   -0.041     0.042  -0.124 0.042   0.333
##              Penalized IVW   -0.010     0.016  -0.042 0.022   0.543
##                 Robust IVW   -0.007     0.018  -0.042 0.028   0.690
##       Penalized robust IVW   -0.012     0.016  -0.043 0.020   0.469
##                                                                    
##                   MR-Egger   -0.113     0.100  -0.309 0.084   0.260
##                (intercept)    0.006     0.007  -0.009 0.021   0.428
##         Penalized MR-Egger   -0.048     0.038  -0.122 0.027   0.213
##                (intercept)    0.003     0.003  -0.003 0.009   0.297
##            Robust MR-Egger   -0.020     0.039  -0.097 0.056   0.603
##                (intercept)    0.001     0.003  -0.004 0.006   0.690
##  Penalized robust MR-Egger   -0.045     0.032  -0.108 0.017   0.155
##                (intercept)    0.003     0.002  -0.002 0.007   0.274
mr_plot(MRInputObject_prostate,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)

MR-PRESSO prostate

#MR-PRESSO prostate

input = data.frame(betaY=diabetes_p$Effect_p,seY=diabetes_p$StdErr_p,beta_diabetes=diabetes_p$Effect,se_diabetes=diabetes_p$StdErr, row.names= NULL)

mrpresso_out = mr_presso(BetaOutcome = "betaY", BetaExposure = "beta_diabetes", SdOutcome = "seY", SdExposure="se_diabetes", data=input, OUTLIERtest = TRUE, DISTORTIONtest = TRUE, NbDistribution = 1000)
## Warning in mr_presso(BetaOutcome = "betaY", BetaExposure =
## "beta_diabetes", : Outlier test unstable. The significance threshold of
## 0.05 for the outlier test is not achievable with only 1000 to compute
## the null distribution. The current precision is <0.088. Increase
## NbDistribution.
mrpresso_out$`Main MR results`
##        Exposure       MR Analysis Causal Estimate         Sd     T-stat
## 1 beta_diabetes               Raw     -0.04085910 0.04223478 -0.9674276
## 2 beta_diabetes Outlier-corrected     -0.01429709 0.01722820 -0.8298654
##     P-value
## 1 0.3360120
## 2 0.4089954

MR-PRESSO Distortion test

mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Distortion Coefficient`
## beta_diabetes 
##     -185.7861
mrpresso_out$`MR-PRESSO results`$`Distortion Test`$Pvalue
## [1] 0.044

MR-PRESSO Outliers detected

out_indices = mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Outliers Indices`
outliers=diabetes_p$SNP[out_indices]


ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%outliers,yes='red',no='blue'))) + 
  geom_point(size=1.5) +
  labs(title='',x = "Genetic association with T2D", y = "Genetic association with prostate cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "outliers"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 10), axis.title.x = element_text(size = 12),
        axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 12),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

MR stratified analysis

BMI

BMI=read.csv('C:/Users/marta/Desktop/IC/Thesis/BMI.csv', stringsAsFactors = FALSE, header=FALSE)
index<-match(diabetes_p$SNP,BMI$V2)
data_p_bmi<-na.omit(diabetes_p[index,])

Creating scatter plot

ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%BMI$V2,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "BMI SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

MRInputObject_prostate_bmi <- mr_input(bx = data_p_bmi$Effect,
                          bxse = data_p_bmi$StdErr,
                          by = data_p_bmi$Effect_p,
                          byse = data_p_bmi$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=data_p_bmi$SNP)

MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_bmi, method = "all")
MRAllObject_all_prostate 
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median   -0.019     0.034  -0.085 0.047   0.569
##            Weighted median   -0.036     0.032  -0.099 0.026   0.257
##  Penalized weighted median   -0.048     0.032  -0.111 0.015   0.135
##                                                                    
##                        IVW   -0.016     0.027  -0.070 0.038   0.561
##              Penalized IVW   -0.032     0.024  -0.078 0.015   0.179
##                 Robust IVW   -0.026     0.026  -0.077 0.025   0.322
##       Penalized robust IVW   -0.034     0.023  -0.079 0.011   0.139
##                                                                    
##                   MR-Egger   -0.006     0.073  -0.150 0.138   0.934
##                (intercept)   -0.001     0.006  -0.013 0.011   0.884
##         Penalized MR-Egger   -0.050     0.064  -0.175 0.075   0.433
##                (intercept)    0.002     0.005  -0.009 0.012   0.759
##            Robust MR-Egger   -0.027     0.053  -0.132 0.078   0.611
##                (intercept)    0.000     0.004  -0.008 0.008   0.961
##  Penalized robust MR-Egger   -0.049     0.042  -0.132 0.034   0.244
##                (intercept)    0.001     0.004  -0.006 0.009   0.701
mr_plot(MRInputObject_prostate_bmi,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)

Insulin

index<-match(diabetes_p$SNP,Insulin$V2)
data_p_insulin<-na.omit(diabetes_p[index,])
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%Insulin$V2,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "Insulin SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

MRInputObject_prostate_insulin <- mr_input(bx = data_p_insulin$Effect,
                          bxse = data_p_insulin$StdErr,
                          by = data_p_insulin$Effect_p,
                          byse = data_p_insulin$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=data_p_insulin$SNP)

MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_insulin, method = "all")
## Warning in lmrob.fit(x, y, control, init = init): M-step did NOT converge.
## Returning unconverged SM-estimate
MRAllObject_all_prostate 
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median    0.084     0.096  -0.103 0.272   0.378
##            Weighted median    0.088     0.078  -0.066 0.242   0.262
##  Penalized weighted median    0.088     0.083  -0.074 0.250   0.288
##                                                                    
##                        IVW    0.078     0.081  -0.080 0.236   0.333
##              Penalized IVW    0.098     0.060  -0.019 0.215   0.100
##                 Robust IVW    0.091     0.232  -0.364 0.547   0.694
##       Penalized robust IVW    0.090     0.064  -0.035 0.216   0.159
##                                                                    
##                   MR-Egger    0.069     0.162  -0.249 0.388   0.669
##                (intercept)    0.001     0.011  -0.021 0.022   0.949
##         Penalized MR-Egger    0.068     0.122  -0.170 0.306   0.577
##                (intercept)    0.001     0.007  -0.013 0.014   0.914
##            Robust MR-Egger    0.048        NA      NA    NA      NA
##                (intercept)    0.002        NA      NA    NA      NA
##  Penalized robust MR-Egger    0.066     0.126  -0.180 0.312   0.599
##                (intercept)    0.001     0.005  -0.009 0.011   0.870
mr_plot(MRInputObject_prostate_insulin,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)

Glucose

index<-match(diabetes_p$SNP,Glucose$V2)
data_p_glucose<-na.omit(diabetes_p[index,])
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%Glucose$V2,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "Glucose SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

MRInputObject_prostate_glucose <- mr_input(bx = data_p_glucose$Effect,
                          bxse = data_p_glucose$StdErr,
                          by = data_p_glucose$Effect_p,
                          byse = data_p_glucose$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=data_p_glucose$SNP)

MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_glucose, method = "all")
MRAllObject_all_prostate 
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median    0.060     0.107  -0.148 0.269   0.571
##            Weighted median    0.065     0.080  -0.093 0.222   0.420
##  Penalized weighted median    0.027     0.077  -0.125 0.179   0.727
##                                                                    
##                        IVW    0.135     0.064   0.009 0.261   0.036
##              Penalized IVW    0.067     0.060  -0.051 0.185   0.266
##                 Robust IVW    0.036     0.063  -0.087 0.160   0.562
##       Penalized robust IVW    0.055     0.062  -0.066 0.176   0.376
##                                                                    
##                   MR-Egger    0.238     0.110   0.023 0.453   0.030
##                (intercept)   -0.007     0.006  -0.020 0.005   0.252
##         Penalized MR-Egger    0.238     0.110   0.023 0.453   0.030
##                (intercept)   -0.007     0.006  -0.020 0.005   0.252
##            Robust MR-Egger    0.167     0.411  -0.640 0.973   0.685
##                (intercept)   -0.005     0.014  -0.033 0.023   0.732
##  Penalized robust MR-Egger    0.167     0.411  -0.640 0.973   0.685
##                (intercept)   -0.005     0.014  -0.033 0.023   0.732
mr_plot(MRInputObject_prostate_glucose,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)

Lipids

index<-match(diabetes_p$SNP,Lipids$V2)
data_p_lipid<-na.omit(diabetes_p[index,])
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%Lipids$V2,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "Lipid SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

MRInputObject_prostate_lipid <- mr_input(bx = data_p_lipid$Effect,
                          bxse = data_p_lipid$StdErr,
                          by = data_p_lipid$Effect_p,
                          byse = data_p_lipid$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=data_p_lipid$SNP)

MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_lipid, method = "all")
MRAllObject_all_prostate 
##                     Method Estimate Std Error 95% CI         P-value
##              Simple median    0.075     0.103  -0.127  0.277   0.465
##            Weighted median    0.083     0.093  -0.099  0.266   0.372
##  Penalized weighted median    0.053     0.090  -0.124  0.229   0.558
##                                                                     
##                        IVW    0.154     0.074   0.008  0.300   0.038
##              Penalized IVW    0.090     0.067  -0.042  0.222   0.182
##                 Robust IVW    0.114     0.173  -0.225  0.453   0.508
##       Penalized robust IVW    0.081     0.068  -0.053  0.215   0.238
##                                                                     
##                   MR-Egger    0.417     0.120   0.182  0.653   0.001
##                (intercept)   -0.019     0.008  -0.034 -0.004   0.012
##         Penalized MR-Egger    0.417     0.120   0.182  0.653   0.001
##                (intercept)   -0.019     0.008  -0.034 -0.004   0.012
##            Robust MR-Egger    0.420     0.122   0.181  0.660   0.001
##                (intercept)   -0.019     0.007  -0.033 -0.006   0.006
##  Penalized robust MR-Egger    0.420     0.122   0.181  0.660   0.001
##                (intercept)   -0.019     0.007  -0.033 -0.006   0.006
mr_plot(MRInputObject_prostate_lipid,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)